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Abstract 

The split-operator pseudo-spectral method based on the fast Fourier transform (SO-FFT) is 
a fast and accurate method for the numerical solution of the time-dependent Schrodinger-like 
equations (TDSE). As well as other grid-based approaches, SO-FFT encounters a problem of the 
unphysical reflection of the wave function from the grid boundaries. Exterior complex scaling (ECS) 
is an effective method widely applied for the suppression of the unphysical reflection. However, 
SO-EFT and ECS have not been used together heretofore because of the kinetic energy operator 
coordinate dependence that appears in ECS applying. We propose an approach for the combining 
the ECS with SO-FET for the purpose of the solution of TDSE with outgoing-wave boundary 
conditions. Also, we propose an effective ECS-friendly EFT-based preconditioner for the solution 
of the stationary Schrodinger equation by means of the preconditioned conjugate gradients method. 

PACS numbers: 02.60.Cb, 02.70.Hm, 02.70.Bf, 02.60.Dc, 32.80.-t, 33.80.Gj 
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I. INTRODUCTION 


A numerical approach based upon the split-operator method and the fast Fourier trans¬ 
form (SO-FFT) is widely used for the solution of equations of the form of time-dependent 
Schrodinger equation (TDSE). Such equations commonly appear in the consideration of the 
wide range of problems of a quantum mechanics and an optics. So far SO-FFT approach 
has been applied for the aim of the solution of one- and two-dimensional TDSE in the 


Cartesian coordinates 


fl. 


three-dimensional TDSE in the spherical coordinates 2|, l3| , the 


time-dependent Kohn-Sham equations the Gross-Pitaevskii equation j^, and molecular 
dynamics joj as well. The method is non-effective in use for singular potentials ^ because 
of both the inaccuracy of the Fourier transform approximation for the wavefunction with 
discontinuous derivative and the split operator method conditional stability. But the kind 
of problems characterized both by a smooth potential and a spectrum strictly bounded from 
above imply the Fourier transform to converge more rapidly jb] by contrast to the hnite 
difference schemes, the expansions over B-splines and hnite elements, and any other polyno¬ 
mial approximation based approaches. As compared to other time propagation techniques 
commonly used along with FFT, namely the leap-frog method and short iterative Lanczos 
method, the split-operator method appears to demand the less amount of memory for the 


0 . 


service arrays storage PJ. The latter is particularly important for the many-dimensional 
problems. 

As a grid based approach SO-FFT may provide the solution only on a hnite space region. 
If the wavefunction leaking beyond its boundaries becomes apparent during the evolution, 
then one encounters the unphysical rehection from the grid boundaries. Upon using SO- 


FFT the unphysical rehection is commonly suppressed 


2| with the help of the adding of 


the complex absorbing potential (CAP) js]. However the CAP method exhibits a number 
of shortages. That are the complex potential absorbs best the small energy wavefunction 
components while the grid boundary is most fast reached by the high energy components. 


and also it distorts the wavefunction behavior in the CAP absent region 8l4ll|. 


The CAP method disadvantages are absent in the exterior complex scaling (ECS) ap¬ 
proach 12|]. For the short iterative Arnoldi-Lanczos method the way for combining of the 
FFT and ECS has been demonstrated in |^. However the application of SO-FFT along 
with ECS is impeded by the explicit coordinate dependence of the scaled kinetic energy 
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operator [^. In the present work a techniqne for the snperposing SO-FFT and ECS is snp- 
posed. The technique is based upon the splitting of the scaled kinetic energy operator to the 
spatially homogeneous part (for which a step may be performed by means of FFT) and the 
absorbing term (which can be handled with the help of the hnite difference). The approach 
under proposal has been tested through a calculation of continuum spectrum amplitudes for 
the model system by t&E-SURFF method j^. For the purpose of the iterative solution of 
the stationary Schrodinger equation (SSE) (which is necessary for the accomplishing of the 
E-SURFF part of the t&E-SURFF scheme) we suggest the FFT-based preconditioner. 

The paper is organized as follows. In the Section [ITl the proposed method for the SO-FFT 
and ECS combining is exposed. The Section [ITIl describes the using of the preconditioned 
conjugated gradient method (PCG) with the aid of performing the SSE iterative solution, 
and also develops the advanced FFT-based preconditioner designed for the using with ECS. 
The section HVl presents the SO-FFT testing, including the proposition of the advanced ECS 
contour. 

Below we use the atomic units of measurement. 

II. SOLUTION OF THE COMPLEX SCALED TDSE BY MEANS OF SO-FFT 


Let us consider TDSE 


5'0(x, t) 

where the kinetic energy operator 


K + U{x^ t) ip{x^ t), 


K = 


( 1 ) 


2 dx^ 


( 2 ) 


The formal solution of the TDSE is 


^/J{x,t + T) = e 


( 3 ) 


The solution can be implemented by using constant-step discrete variable representation 
(DVR) l4] over a space variable. Let us set up in the region x G [a, b] the uniform grid 


Xi = a + h{i — 1/2), i = 1,..., A, 


( 4 ) 
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where iV = 2*^ is the nodes number (z/ being an positive integer) and h = {b — a)/N is the 
grid step. We can expand 'ijj{x,t) over the basis of Lagrange functions 

N 

^i{x) = (5) 

n=l 

Here 


^n{x) = 


exp{iknx) 
\/b — a 


are the Fourier basis functions, where the momentum kn is 


k 


n 


ne|l.iV/2]i 

+ ne|iv/2 + i,iV]. 


( 6 ) 


It is easy to check that Lagrange functions have property ii{xj) = VhSij^ and coefficients of 
'ip{x,t) expansion over Lagrange functions are = Vh'ijj{xi,t). 

The evolution of the function V’(^) values vector at the grid nodes may be evaluated with 
the help of the split operator method along with FFT as follows: 




(7) 


Here F and F’ -1 = Ft are the orthogonal matrixes of the direct and the inverse discrete 
Fourier transforms (FT) respectively. Their elements are defined by the expression 

Fni ( 8 ) 


Due to the transition to the momentum representation through the direct FT and to the 
DVR through the inverse FT the matrixes of the potential U and the kinetic K energy are 


diagonal, that is Uij = U{xi,t + T/2)6ij and K^m = ^^ij- Thus the action of the propagators 
and reduces to the simple wavefunction multiplication by the phase factors. The 
multiplication by a fully populated matrix F generally requires computational operations. 
But it might be performed by using FFT that requires the number of operations of the order 
of N log 2 N. 

In order to provide the outgoing wave boundary conditions one may use the ECS method 
consisting in the change of the real space variable to the complex one 
way: 


12l | in the following 


dt 


1 


2 dz"^ 


+ U{z,t) 




(9) 
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Here the integration complex contour is conveniently described 12|, ll6| by expression 


z{x) = / q{x)dx. 

Jo 


( 10 ) 


Let us introduce the new function 


ij^x) = ^/q{x)'ijj{z{x)) 


( 11 ) 


and express the differentiation over the complex z through the differentiation over the real 
X as 


diplx, t) 
dt 


did 


2^/q{x) dx q{x) dx 


+ U{z{x),t) ) ^{x,t). 


( 12 ) 


It could be seen from here that upon the complex scaling performing the direct applying of 
FFT becomes impossible because the appearing of the explicit x dependence in the kinetic 
energy operator 


K, = 


did 


2^/q{x) dx q{x) dx 

Hence the functions (E]) would not be its eigenfunctions anymore. 

But we may divide the Kg operator into the two terms as Kg = K + C^ where 


(13) 


c = k,- k. ( 14 ) 

By construction C = 0 in the domain where q{x) = 1, i.e. outside the ECS region. For 
the purpose of the operator K diagonalization one may use FFT, while the operator C 
(that actually provides the wavefunction absorbing in the scaling region) approximation can 
be performed through the hnite-difference scheme since high accuracy calculation of the 
wavefunction is not necessarily needed in the ECS region. 

We shall introduce the matrix C as the difference 


C = K, - Kfd, 


(15) 


where and Kpo are respectively the operators Kg and K approximations implemented 
by the symmetric three-point difference scheme as 

.-V2 


[K.v>l. = -F' 


2 h 


- Qi Qi 


qi+i/2h 


qi-i/2h 


[KpoV^li — 


1 V'i+i - 2i)i -f 

'2 K 


(16) 

(17) 
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This leads to a complex symmetric tridiagonal matrix C. The matrix C symmetry is gen¬ 
erally not crucial in a TDSE solution, but this property is needed for the providing of PCG 
convergence in a SSE solution process (that is considered in the next Section). 

Thus the numerical scheme takes the following form: 

+ r) = (18) 


where the multiplication by the exponential function of the tridiagonal matrix C may 
be performed with the help of the Cranck-Nicholson (CN) scheme as 

/ • \ -1 / • \ 


^-iCr 


-ll + yC 




(19) 


CN requires the number of operations of the order of N. 


III. SOLUTION OF THE STATIONARY SCHRODINGER EQUATION USING 
EFT PRECONDITIONER 


The i mp 


SURFF 
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ementation of the E-SURFF approach and its generalization known as t&E- 
both require the solution of the SSE with the right-hand side of the form 


ks U{z{x),tr^ 


-E 


^e{x) = f{x) 


( 20 ) 


Upon using the approximations described in the previous Section the operator in the left- 
hand side transforms into the matrix 


A = F-^KF + C + U - FI. (21) 

Here I is the identity matrix, Jjj = 5ij. The matrix A is a fully populated, so the direct 
solution of the equation 

At/? = f (22) 

by the Gaussian elimination method requires a number of operations of the order of N^. 
However multiplication of a vector by the matrix A may be performed through N log 2 N 
operations since the main overhead is caused by the multiplication by F~^KF. The multi¬ 
plication by the tridiagonal matrix C as well as by the diagonal matrixes requires a number 
of operations of the order of N. 
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Since the matrix A is complex symmetric the system of equations (|2^ might be solved 
by means of an iteration method namely the preconditioned conjugated gradients (PCG) 


method 


18 


19[ |. On every PCG iteration one performs the multiplication of a vector by 
the matrix P^^A, where P is a preconditioner matrix that satishes the equation P^^A I 
yet at the same time leads to a small operations number required for the equation Py = r 
solving. 

In our case the natural choice is a preconditioner based on FFT in the following way 


p_i ^ F-1(K- 


(23) 


where 

[(K - ET)-\^ =(^-E^ 5nm. (24) 

Upon utilizing the CAP approach (which is covered in details in the next Section) the simple 
preconditioner fl2^ provides the fast PCG convergence. ECS test calculations showed that 
upon the condition |g(a;)| ~ 1 PCG converges by rather small iterations number. But for 
the TFECS approach (see the next Section) q{x) —)■ oo near the grid boundaries that leads 
to the extremely slow convergence. 

For this reason especially to deal with ECS we have developed the advanced FFT-based 
preconditioner as follows. Let us split the spectrum into the bands each containing the set 
of the momenta /c„, n = ri;,..., u; + A; — 1. Then let us write down an eigenvalue problem 
for the square matrix acting on the vector components in the one band only as follows 


Here Xun 


0 at n ^ [n ^,ni + Ni — 1], the solution number u = ni,... ,ni + Ni — 1, and 


(25) 


1K,| 


nm 


^hnkmQn—ni+l,m—ni+li 777. G T A; l], 

0, n,m ^ [ni,ni + A) - 1]. 


(26) 


In turn, here 

Qnm= [ (Pn{x)[q{x)]~^(pm{x)dx. (27) 

J a 

Since q{—x) = q{x), true is Qmn = Qnm, therefore the matrix is complex symmetric. 
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The block diagonal matrix K = be considered as the matrix of an operator 


Kr. 


Id 1 d 

2 dx [q{x)Y dx' 


(28) 


with artificially nnllified off-blocks elements. Using Ka instead of Kg here yields the following 
advantage. This choice allows to calculate Qmn elements only once for 1 < n, m < Ni, 
then matrix K/ elements for any n,m ^ obtained from the Qmn through the 

simple multiplication by knkm/^, according to fl26|) . The distortion of the eigenvalues and 
eigenvectors caused by the distinction of Ka from Kg are unessential by compare to the 
distortions because of the off-block elements nullifying. 

Further let us construct the matrix X by using the equation solutions as the rows 
in the following way: 


Xvm = [xAL- (29) 

Since the matrix K is complex symmetric then the matrix X should be complex orthogonal 
that is X-i = X^. 

The preconditioner we are proposing here is 


p-i = F-1X-1(A - Eiy^XF, 


(30) 


where A is a diagonal matrix with elements where Xy is an eigenvalue from 

fl25|l . The matrix X is block-diagonal, therefore multiplication by it requires the operations 
number of the order of NiN. Thus upon the condition A) < log 2 N the solution of the 
preconditioner equation under proposal would not be essentially more time demanding than 
the one for the usual FFT preconditioner. At that the PCG convergence speeds up essentially 
due to the proximity of eigenvalues and eigenvectors of P for the ones of the matrix A. 

The Fig. [1] displays several typical preconditioner basis functions in the coordinate rep¬ 
resentation Xv{xi) = [F“^Xj,]j for the case Ni = 10. These functions break down into the 
two kinds: the ones localized in the ECS range (upper row of Fig. [1]) and the ones dif¬ 
fused over the non-scaled region (lower row of Fig. [T]). The former typically have a large 
negative imaginary part of its eigenvalue. Since by construction the Xu{^) have momentum 
uncertainty Ak < -^Ni then the minimal region of its possible localization has the size 
Ax > {b — a)/Ni due to the uncertainty principle. Thus one can draw a condition on the 
minimal Ni needed for at least one function to be localized in the scaling region as follows: 



A = 3.1876 - 2.35861 


A = 10.618-2.28581 



- 0.1 I i . ; 

-50 -25 0 25 50 

x(a.u.) 





FIG. 1: Typical preconditioner basis functions with corresponding eigenvalues. Solid lines — 
absolute value, dotted lines — imaginary part. 


Ni > {b — a)/2Axcs, where 2Aa;cs is a joint width of the two scaling regions located at the 
grid edges (as we use the basis functions (E]) that are periodic, so these two regions may be 
considered as stitched together through the grid external boundaries). 

It is not possible to split the spectrum into the bands of equal width Ni = const in a 
kn symmetric way, while ki = 0 and kN/ 2+1 has undehned sign, so the including of the 
corresponding components in any band leads to the disparity of the bands associated to the 
positive and negative kn- Therefore we used the following splitting: two bands of the width 
Ni = 1 each containing only one component, namely ki and kjy/ 2 + 1 , and the set of constant 
width Ni = IVbw = 10 bands, and hnally the two mutually symmetric bands of the width 
Nbw a Ni < 2Nbw ~ 1 which include the rest kn values on either side from kj^/ 2 +i- 

It should be noted that the FFT preconditioner (in both proposed variations) is easily 
generalized to the many-dimensional case. In j^, ^ the solution of the many-dimensional 
SSE has been performed by means of the preconditioner based on the direct solution of 
the sparse block-tridiagonal matrix system. Authors of 2l|] have used the preconditioner 
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based on the orthogonal transformation to the representation of the one-dimensional problem 
eigenvectors set. Either of these two preconditioners utilizing implies carrying the number 
of operations of the order of (here N being the full number of the two-dimensional grid 


points). FFT-based preconditioner needs the number of operations o 


only for same case. It is inferior to the multigrid preconditioner 


22 


he order of N log 2 N 


23| . requiring as few 


as N operations. However the FFT approach over-performs the multigrid approach when 
used for the problems for which the FFT a.pproach exceeds the hnite-difference schemes in 
the space approximation convergence rate 241. 


IV. TESTING OF THE PROPOSED METHODS 


As a model example we have chosen the Gaussian packet 

„2 


'0(a:, 0) = TT 


= ^- 1/4 


X 


exp ( —-+ ivqX 


moving with the speed uq = 2 in the free space. The Gaussian packet spectrum is 

{k-v^r 


^teorik') ^ 


= ^-1/4 


exp 


(31) 


(32) 


So we can test the methods for the unphysical reflection suppression by comparing the 
theoretical spectrum to the one obtained by means of t&E-SURFF. 

For all the examples listed below we used the same grid parameters as follows. The grid 
boundaries were supposed to be 6 = —a = r^ax = 51.2, the EGS region was |x| > Vg = 40. 
The time step was taken to be equal to r = h being the space grid step. The evolution 
has been performed up to the time tmax = 20 (unless otherwised noted). At this t^ax the 
packet goes into the scaling region exactly by half (see the Figj2]). Such a choice allows 
to estimate simultaneously the accuracy of the both t&E-SURFF components, namely t- 
SURFF (which implies an amplitude extraction from the probability amplitude flux passed 
through the scaling region boundary) and E-SURFF (that extracts an amplitude from the 
rest function), since at such t both components contribute equally to an amplitude under 
evaluation. 

The EGS is usually introduced by the sharp contour rotation to the complex plane as 
follows: 

1, \x\ < Tg] 


q{x) = 


r,ida 


\x\ > Tg, 
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x(a.u.) 

FIG. 2: Wavefunction for t = 20 calculated by using SO-FFT with SECS. Solid line — absolute 
value, dashed line — imaginary part. 


where Og is the angle of the rotation to the complex plane, and Vg is the distance from the 
origin to the rotation onset points. However npon snch a choice C has a discontinnity at 
the scaling region bonndary that leads to the loss of the fnnction smoothness at the 
points I a: I = Vg. For the non-smooth fnnctions the Fonrier expansion converges extremely 
slowly. It manifests itself in the strong wavefnnction reflection from the scaling region. 



x(a.u.) 


x(a.u.) 


(a) 


(b) 


FIG. 3: The ECS contours under consideration: left — SECS, right — TEECS. Solid lines — 
imaginary part, dotted lines — real part. 


Therefore we have used the smooth ECS contour (smooth ECS, SECS) for which (see 
Fig. |3Ka)) 


q{\x\ > Vg) = exp 


i'll |a;| — Vg 

4 Aa:cs 


(33) 
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where the scaling region width Aa:cs = ^max — fs- For such contour, a maximal rotation 
angle is 7r/4 at a; = r^ax- 


CAP (Ecap=2) 


SECS 




TRAP {k„,i„=0.25) 


TFECS 




FIG. 4: The amplitude A{k) (|A| — solid lines, AT — dashed lines): a) CAP; b) SECS; c) TFAP; 
d) TFECS. 


The result of the test amplitude evaluation by means of SECS for the points number 
N = 1024 (corresponding to the grid step h = 0.1) is shown in the FigjU^b). It is seen that 
the calculated curve differs the most from the Gaussian shape at the small momentum k 
values. Under the ECS applying the boundary reflection coefficient has the form 

Recs{E) = exp [-2k^z{rmax)] ■ (34) 

Its tending to 1 at fc —)■ 0 results in the absence of the unphysical reflection suppression for 
the small fc, therefore the situation observed was expected. 

In order to examine the correctness of the phase evaluation we represent in the FigJUalso 
the amplitude imaginary part that should be equal to zero due to fl32|) for the model system 
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under consideration. 

The authors who utilize SO-FFT commonly use the complex absorbing potential (CAP) 
for suppression of the unphysical reflection. It makes sense to compare the CAP and ECS 


efficiencies. As it is shown in 


13|, for the hxed kinetic energy value E = Ecap the ECS 


solution is equivalent to the solution with the complex absorbing potential of the form 

Ua{x) = Ecap[1 - q{xf]- (35) 

As a hrst CAP version we used the potential given by the formula fl3^ for q{x) from 
fl33|) at Ecap = 2 a.u., that is the energy corresponding to the amplitude maximum. Such 
CAP is close to the commonly used linear CAP. The result is presented in the FigJU^a). 
One can easily see the essential curve distortions exceeding those for SECS (Fig. 111(b)). By 
virtue of the quasiclassical approximation it is not difficult to derive the formula for the grid 
boundary reflection under the CAP action: 


Rcap{E) ~ exp 


'^Ua{x)dx 


(36) 


That is the maximum reflection is to occur at the large k values. As the greatest distortions 
in the Fig. Hla) are observed at small k values, they are expected to be caused not by the 
grid boundary reflection, but by the reflection from the CAP itself. 


In the work 


ll| an advanced CAP version is proposed, namely so-called transmission-free 


absorbing potential (TFAP) having the form 


Ua{\x\ > Vs) = -i 


^2 . 

• ^mir 


\x\ — r. 


Ax 


cs 


(37) 


where 


y{s) = 


as — bs"^ + 


(38) 


(1 - s)2 (1 + S)2J ’ 

and parameters are a = — 16, 6 = — 17, c ~ 2.622057. This potential was obtained by 

the authors of [n| from the condition of minimization of the wavefunction deviation from 
the quasiclassical outgoing wave hence the minimization of the reflection from the potential 
itself. The parameter /cmin in 1137)) has a meaning of the smallest momentum needed for the 
reflection from the potential to be small. Upon kmin value increasing the reflection from 
the potential grows, and upon its decreasing one encounters the reflection from the grid 
boundary since the absorbing potential weakens at small kmin- The Fig. m displays the 
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results for = 0.25 a.u. It is apparent that the results are much more accurate than 
those obtained by the linear CAP. 

By analogy with TFAP we propose transmission-free ECS (TFECS) approach making 
use of the following contour: 


where 


q{\x\ > Ts) 



(39) 


y{s) = a 


1 

(T^ 


1 

(1 + S )2 



(40) 


That is we chose q{x) meaning that is equivalent to the absorbing potential of the form 
Ua{x) = —iy [|a; — TsI/Axcs] due to fl35ll . The function y{s) asymptotic behavior at s —)■ 1 
(x —)■ Tmax) E analogous to that for y{s) from (l38|l . however y{s) ~ at s —?■ 0 (x ^ Vg) 
while y{s) ~ s. The reasons for such a choice will be described below. At such q{x) the 
complex coordinate z(x) tends to the inhnity logarithmically at x ^ rmax- But by virtue of 
the grid Xj shift by h/2 the value z{xn) = —z{xi) ~ \/aiAxcs ln(2Axcs/h) becomes hnite 
(see Fig. [3](b)). This is the reason for the choice of Eq.(jl]) for grid dehnition. The Fig. 111(d) 
represents the TFECS results for the parameter a = 0.5. It is seen that TFECS yields much 
more accurate results than those obtained by the rest methods under consideration. 

Upon the Az(xjv) value substitution in the Eq. fl5T)l . one may evaluate the coefficient of 
the grid boundary reflection under TFECS using in the following way: 


Rtfecs{E) ~ 


h 


) \/2akAxQQ 


(41) 


For some physical problems the calculation accuracy is crucial at large k values. Therefore 
we shall consider this issue. The Figj5] shows the results obtained by TFAP, SECS, and 
TFECS in the logarithmic scale. Since the wavepacket moves to the right, the reflection from 
the TFAP/ECS region produces a packet moving to the left that results in the amplitude 
distortion at A: < 0. It is apparent that even the simple SECS suppresses the unphysical 
reflection much better than TFAP does. Nevertheless at fc < —2.5 one can easily discern 
the error of the order of 10“^ entailed by the reflection. As the “exact” complex scaling 
should not cause any reflection at large \k\, so the main reason for the results to deviate 
from the theoretical form at large \k\ appears to be the inaccuracy of the hnite-difference 
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FIG. 5: The absolute value of the amplitude A{k): TFAP (dashed line), SECS (dotted line), and 
TFECS (solid line). 

approximation of the absorbing operator C that leads to a small parasite potential from 
which the reflection does occur. The greatest error arises in the points of the function q{x) 
higher derivatives discontinuity, namely the point of the rotation to the complex plane, 

X = Vs. 

As under the TFECS using the function q{x) has the derivative discontinuity of the 
fourth order only near this point, so TFECS applying results in the reflection from the 
error-caused parasite potential yielding the amplitude error of the order of 10“® only (that 
corresponds to the cross section error 10 ^^). Thus upon the TFECS utilizing ^z{rmax) 
provides the essential unphysical boundary reflection suppression, while the smoothness 
of the complex plane rotation ensures the smallness of the error caused by the absorbing 
operator C approximation by the hnite difference scheme. 

Further let us examine more thoroughly the solution convergence to the exact one with 
the grid step h decreasing, hereinafter under the TFECS using. Since for the aim of ECS 
introduction to SO-FFT we use the three-point finite-difference scheme, it makes sense to 
compare our approach accuracy to that obtained by the three-point hnite-difference scheme 
(fT6|) for the full Kg. We shall use “pure” E-SURFF (that is the amplitude extraction from 
the wavefunction at f = 0) in order to separate the spatial approximation error effect. The 
Figl6](a) presents the error 5{k) = \A{k) — Aeor(^)| of the amplitude obtained by the finite- 
difference scheme, for the spatial step values h = 0.2, 0.1, 0.05. In order to demonstrate 
the relative value of the error, we display on the same hgure the theoretical curve |y4teor(^)|- 
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ESURFF-FD 


ESURFF-FFT 




FIG. 6: The error S = \A{k) — ^teor(fc)| calculated by using E-SURFF (from the wavefunction for 
t = 0) for the different values of the grid points number: N = 512 (dotted line), N = 1024 (dashed 
line), N = 2048 (solid line). The dependence |74teor(fc)| is also shown (thin solid line). Left — FD, 
right — FFT. 

One can easily see the error largeness. As expected for the second-order scheme, the error 
decreases by the factor of 4 upon the step twofold decreasing. 

The Figl6](b) represents the result obtained with the help of FFT. It is apparent that the 
error is many orders smaller than that given by FD (except for the one for the smallest values 
k). The main error is entailed by the parasite reflection arising due to the inaccuracy of the 
hnite-difference C approximation. Since we use the hnite-difference scheme of the second 
order, the coefficient of the parasite reflection should depend on h quadratically, exactly as 
it is observed in the Figl6](b). Though formally the convergence rate is the same as upon 
the using of the hnite-difference scheme for the full Kg, the absolute value of the error is 
much smaller. 

The peak of the error at /c = 0 arises due to the rehection from the grid boundary. There¬ 
fore, this peak is hardly reduced with decreasing of h. It may be suppressed signihcantly by 
the scaling region Axes increasing only, as it follows from Eg. (1411) . 

The Fig. [7] shows 6{k) for SO-FFT at tmax = 20 a.u. The error increasing as compared 
to that of E-SURFF (Figj6](b)) is mainly caused by the error of the approximate integration 
rule for the time integration of the probability amplitude hux In the Fig. [7](a) one 
can see the error appearing upon the using of Simpson’s rule for the time integration of the 
probability amplitude flux. At fc > 5 a.u. the integration error dominates the other error 
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SO-FFT (Simp.) 


SO-FFT (Filon-Simp.) 




FIG. 7: The error 5 = \A{k) — Ateor(fc)| calculated by using tE-SURFF (from wavefunction for 
t = 20) for the different values of the grid points number: N = 512 (dotted line), N = 1024 (dashed 
line), N = 2048 (solid line), with the time integration performed by Simpson’s rule (left) and with 
Filon-Simpson’s rule (right). The dependence |^teor(^)| is also shown (thin solid line). 

sources. Since we used the time step t = h? whereas the Simpson’s rule error is 0(r'^), so 
the integration error depends on the grid step h as 0{h^), as is observed in the hgure. The 
reason of the error growing upon the \k\ increasing is that the probability amplitude flux at 
large energy E is the rapidly oscillating function on t (that is the product of exp{iEt) and 
quite slowly changing functional of the wavefunction). The Figl7](a) demonstrates the error 
arising upon the using of Filon-Simpson’s rule for the time integration of the probability 
amplitude flux that is the formula especially developed for the fast oscillating functions 
integration. It is seen that in contrast to the Simpson’s rule applying results the error does 
not grow under k increasing in the region of the large \k\ values. But the error grows a bit 
at small \k\. 

Finally let us make some remarks upon the convergence of the PCG approach used for 
E-SURFF. During the TFECS using with N = 1024 and preconditioner of the form fl30|l 
with Ni = 10 the average (over all the calculated spectrum points) iterations number was 
Niter ~ 79 until the discrepancy of the equation fl2^ was equal to £ = 10“^^. In contrast, 
during the SECS using with the same grid and preconditioner parameters the average PCG 
iterations number was as small as Nuer ~ 36. The reason for such a situation is that q{x) 
differs from 1 strongly in some points for TFECS, while it does not for the SECS. Thus 
in the case when small E amplitudes are not needed it makes sense to prefer SECS more 
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likely than TFECS. Whereas the calculation accuracy is not crucial one might utilize even 
the TFAP approach as it implies the PCG iterations number being as small as Wter ~ 20 
under using the preconditioner of the form fl30p . 


V. CONCLUSION 

In the present work we are developing the technique aimed to the combining of FFT and 
ECS approaches for the TDSE solution in the context of the important quantum mechanical 
and optical problems consideration. The technique consists in the TDSE solution by means 
of the split-operator method through the splitting of the kinetic energy operator into the 
coordinate-independent and coordinate-dependent parts. The coordinate-dependent part 
can be approximated by the hnite-difference scheme. 

In the present work the proposed technique efficiency is proved by means of the contin¬ 
uum amplitudes calculation through the t&E-SURFF approach [^. The latter implies the 
solution of not only the TDSE, but of the SSE with the right-hand side as well. So we 
have also demonstrated that the solution of the SSE with the right-hand side by means of 
the conjugate gradient method may be performed with the help of the FFT-based precon¬ 
ditioner. 

Further, during the testing we investigated the ECS method efficiency with different 
contours in use, bearing in mind that the Fourier transform converges rapidly for a smooth 
function only. As a result we may conclude that the ECS approach can be accomplished 
by making use of the advanced TFECS contour proposed here as well. It provides the high 
degree of the unphysical reflection suppression for all energy values along with the small 
reflection from the approximated absorbing operator itself. 

Finally, all the above mentioned actions taken resulted in the error decreasing such that 
the time integration error became distinguishable. So we also examined the efficiency of the 
numerical integration algorithm in the t-SURFF implementation. We can make a conclusion 
that the time integration of the probability amplitude flux is preferably to be performed by 
means of the Filon-Simpson’s rule. 

In summary, we have elaborated the highly efficient and fast method which can be helpful 
in solving a number of quantum mechanical and optical problems involving the evaluation 
of the amplitudes of ionization and dissociation of atoms and molecules by external helds. 
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the laser beams propagation in the waveguides with leaking and so on. Thereupon, in order 
to apply the complex approach developed in the current article to the real problems, in the 
future we plan to upgrade it for a multidimensional case. 
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